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Molecular noise in gene regulatory networks has two intrinsic components, one part being due to fluctuations 
caused by the birth and death of protein or mRNA molecules which are often present in small numbers 
and the other part arising from gene state switching, a single molecule event. Stochastic dynamics of gene 
regulatory circuits appears to be largely responsible for bifurcations into a set of multi-attractor states that 
encode different cell phenotypes. The interplay of dichotomous single molecule gene noise with the nonlinear 
architecture of genetic networks generates rich and complex phenomena. In this paper we elaborate on 
an approximate framework that leads to simple hybrid multi-scale schemes well suited for the quantitative 
exploration of the steady state properties of large-scale cellular genetic circuits.Through a path sum based 
analysis of trajectory statistics we elucidate the connection of these hybrid schemes to the underlying master 
equation and provide a rigorous justification for using dichotomous noise based models to study genetic 
networks. Numerical simulations of circuit models reveal that the contribution of the genetic noise of single 
molecule origin to the total noise is significant for a wide range of kinetic regimes. 


I. INTRODUCTION 

Molecular noise functions as a driving force for phe¬ 
notypic heterogeneity in cell populations^— . With the 
growth of sophisticated imaging methods and biochem¬ 
ical techniques that probe gene expression at a single 
cell level, the omnipresence of molecular noise in cellu¬ 
lar processes has become widely appreciated-^. Noise in 
biology therefore should no longer be viewed as a nui¬ 
sance with which cells must cope but rather it must be 
viewed as a natural facilitator for adaptation, growth 
and development^^— . Molecular noise in genetic net¬ 
works has at least two source s 9 ’ 10 : one source is the 
probabilistic nature of molecular associations such as the 
binding/unbinding of promoters to operator sites on the 
DNA while the other obvious source is randomness in the 
sequential reactive events such as the stochastic synthe¬ 
sis and degradation of mRNA, protein molecules, which 
generate the whole population of molecules. The di¬ 
chotomous Markov noise (DMN) that arises from a single 
copy of a gene stochastically switching between its ON 
and OFF states, whether through promoter binding or 
through architectural changes of chromatin, is of a quali¬ 
tatively different mathematical character from the other 
sources of molecular noise which may be treated as con- 
trollably small perturbations when the molecular popu¬ 
lations of mRNA and proteins are sufficiently larged. 

Dichotomous gene noise leads to bursty intermittent 
production of mRNAs. The intermittency arises be¬ 
cause mRNA molecules are rapidly synthesized in large 
groups when the gene becomes activated but these bursts 
are then followed by intervals of silence when the gene 
is turned off. ’’Poissonian” stochastic kinetics arises 
when the mRNA’s are synthesized independently from 
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each other at random intervals with uniform probability. 
In contrast, single cell experiments provide ample evi¬ 
dence that the statistics of gene expression is often non- 
Poisonian, with bursty behavior having been observed in 
both prokaryotic*!! - — and eukaryotic*!^— cells. As a rule 
the bursty stochastic kinetics leads to higher overall levels 
of noise. The detailed microscopic origins of dichotomous 
gene noise and the observed transcriptional and transla¬ 
tional bursts still remain under debate— Possible mech¬ 
anisms include local chromatin remodeling, thermally in¬ 
duced nucleosome turnover, covalent modifications and 
genome wide chromatin fiber dynamics— 

Dichotomous noise generally breaks detailed bal¬ 
ance in the circuits and thus creates non-equilibrium 
steady states which can not always be described by 
quasi-equilibrium fluctuation statistics. Dichotomous 
noise driven phenomena include robust phase synchro¬ 
nization^^, stochastic hypersensitivit y 24 ’ 25 , enhanced 
stochastic resonance!!, hysteresis!! and patternin g 28 ’ 29 . 

Brute force simulation of the full master equation for 
genetic networks has already yielded many insight s 30 ' 31 . 
Nevertheless in practice the sheer number of cellular com¬ 
ponents often renders such an explicit approach over¬ 
whelming. Most conventional approximations for treat¬ 
ing stochastic chemical reactions, such as size expan¬ 
sion method s 32 ' 33 are not applicable when there is suffi¬ 
cient dichotomous noise since these approximations usu¬ 
ally rely on the uniform scaling with size of the entire 
population of species to the near deterministic limit—. 
Time averaging can play a role however in reducing the 
influence of the single molecule dichotomous noise. In the 
limit of infinitely fast switching of the gene one can again 
obtain the strictly deterministic result if the molecular 
populations of all the other components are sufficiently 
large. This fast gene switching regime has come to be 
known as the adiabatic limit—. In the adiabatic limit 
there is a very large time-scale separation between gene 
switching at the single molecule level and the biochemi- 
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cal synthesis and degradation reactions in the rest of the 
network which can then again be treated as nearly deter¬ 
ministic in the large system size limit. Thus it is most 
important to develop approximations appropriate to the 
non-adiabatic case. 

Recently several schemes for modeling gene networks 
have been proposed which offer various ways of treating 
the stochastic aspects of gene switching^— . Most of 
these schemes still however rely on there being a suffi¬ 
cient degree of time scale separation between gene and 
protein degrees of freedom, which means that switch¬ 
ing dynamics must usually be assumed to be very fast. 
While such an assumption could be a sound one for many 
gene circuits, the intermediate regime where there is no 
such time scale separation is also of interest especially for 
understanding the behavior of mammalian gene circuits. 
In this contribution we present a framework which lifts 
the requirement for the genes to operate near the adia¬ 
batic limit. This framework leads to a hybrid multi-scale 
scheme for numerically exploring potentially large scale 
genetic networks. We demonstrate the usefulness of our 
scheme by analyzing some well known models of simple 
genetic circuits. 


II. RESULTS 

Schematic discussion of trajectory statistics in gene 
networks. 

In this section we discuss the notion of trajectory prob¬ 
abilities for stochastic gene circuits with dichotomous ele¬ 
ments. This notion motivates hybrid simulation schemes 
based on piecewise deterministic systems of equations. 
We illustrate the ideas using a simple one component 
model of a self regulating gene as the prototype for a gene 
network. Due to its simplicity this idealized problem has 
been the subject of numerous theoretical studiesii^~~. 
The state of the circuit is defined by two dynamical vari¬ 
ables: an operator state variable s = {0,1} describing 
whether the gene is active s = 1 (ON) or inactive s = 0 
(OFF) and the number of proteins n = {0,1,2...oo}. 
Since the gene circuit is stochastic, one is interested 
among other things in knowing the probability p s (n,t) 
of having n proteins and the gene being in state s at 
time t. At this level the most satisfying description of 
such a circuit (if one assumes a well stirred mixture of 
molecules) is given by a discrete space master equation: 


dP ^ = G (P(n - 1, t) - P(n, t))+K({n + 1 )P(n + 1, t) - nP{n , t)) + W„P(n, t) (1) 


where P(n,t) = \po(n,t),pi(n,t)) is the state vector 
specifying the probabilities for each gene/protein state. 
The diagonal birth/death matrices G = diag(go, g\) and 
K = diag(k , k ) contain the rates of synthesis g and degra¬ 
dation k which are assumed to be Poisson processes, 
while W„ is the rate matrix that describes gene state 
switching probabilities as a function of protein numbers. 
In the case when the populations of the species are not 
too small (n 1) such a discrete representation de¬ 
fined with integer population number accuracy(Eq. [TJ 
may no longer be absolutely necessary. When the pop¬ 
ulations of species in the genetic network are sufficiently 
large one may instead treat n as a continuous variable 
and use diffusive dynamics for the evolution of pro¬ 
tein population. This approximation is well justified 
by size expansion arguments. Numerical simulations** 
show that Fokker-Planck approximations can be quan¬ 
titatively adequate even when protein numbers are as 
small as 100. In this limit one may reduce the fully 
discrete master equation description to a master equa¬ 
tion in mixed continuous and discrete variables = 

—VJ + W n P(n,t), where the birth/death terms are re¬ 
placed with J = —A(n) — i d n D{n ) corresponding to an 
effective flux of continuous number of molecules To vi¬ 
sualize trajectories generated by such dynamics consider 
a particular realization of the gene circuit’s history that 
goes through a particular sequence of gene state flips: 


a = {so —> si —>• S 2 —> ■■■ —> Sjv}- These are governed by 
the master equation for the state variable conditioned at 
any time by the species population^/(n,i') = W p s (n,t). 
For the present model this conditional master equation 
reads explicitly: 

kp 0 {n,t)\ = /- k of f k on (n ) \ kp 0 (n,t)\ 
\pi(n,t)J \k 0 ff -k on {n)) \pi{n,t)J 

In between gene flipping events the protein population 
will evolve diffusively v = {no — > n\ —>■ ri 2 —>■ ... —>• njv} 
according to one of the following Fokker-Planck equa¬ 
tions: 

d t p 0 (n,t) = [-d n A 0 (n) + \dlD 0 (n)\p 0 {n,t) 

( 3 ) 

d t pi(n,t) = [ - 5„Ai(n) + i^L>i(n)]pi(n,t) 

where the A s (n) and D s (n) stand for the drift and the dif¬ 
fusion coefficients that correspond to the gene state s at 
any time. For the self regulated gene model these terms 
are A s = g s — kn and D s (n) = g s + kn. In this example 
g s is the protein synthesis rate when the gene is in state 
s while k is the rate coefficient corresponding to a first 
order model of protein degradatio n 1 h 46 . Equations 1-2 
completely define the dynamics of the simplest gene cir¬ 
cuit. The system of equations 1 and 2 is known in math¬ 
ematics as a coupled or switched diffusion process^—. 
It features prominently in the problems of molecular mo¬ 
tors, ratchets and gene networks. The coupled diffusion 
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process is analytically intractable for all but the simplest 
systems. 

Now let us consider a particular realization of the 
stochastic trajectory of the gene state and then ask what 
is the probability given this gene state trajectory of ob¬ 
serving a certain sequence of corresponding protein num¬ 
bers (ni, ri 2 , ..riff) at times (t±, £ 2 , ...ijv)- Most experi¬ 
mental studies measure time dependence of protein con¬ 
centrations rather than follow the gene state changes 
explicitly. First let us consider a system with a fixed 
gene state (e.g let’s say s = 1 at all times). Using 
a convenient bra ket notation the conditional probabil¬ 
ity that protein state rij+i follows rit at an earlier time 
can be expressed asp(nj+i,U+i|nj,U) = (ni+i\U(At)\rii) 
where U(At) is the evolution operator of the Fokker- 
Planck equation corresponding to the particular fixed 
gene state. If we then allow the gene state to change 
in between gene switching events the diffusion of protein 
numbers no n\ —4 ... will be governed by one of 
the Fokker Planck equations (EqO. The evolution op¬ 
erator U s (At) will therefore depend on the gene states 


s at each time t as indicated by its subscript. Thus for 
any particular sequence of gene states a we have propa¬ 
gators p Si (n i+1 ti +1 \ni,ti) = (n i+ i\U Si \ni) defined on the 
intervals between gene switching events. Similarly for 
a particular realization of protein numbers v at times 
ti we can write the probabilities of the gene switching, 
Pm(si+iti+i\si,ti) = (sj + i|G ni |sj). Again this is defined 
on the intervals between protein number diffusion events. 
Suppose we are given set of protein numbers observed 
at regular time intervals At , v = {no, n\, 712 , ...njv} 
which could be obtained by recording a stochastic tra¬ 
jectory from a numerical simulation on a computer or 
by observing the gene expression profile in a single cell 
in the laboratory. Using the assumed Markovian prop¬ 
erty of the underlying molecular stochastic processes we 
can write the probability of observing a trajectory v 
in terms of the protein states irrespective of the gene 
trajectory as the product of transition probabilities be¬ 
tween protein and gene states, which is then summed 
over all possible realizations of the gene state trajectory 
O = {sQj Si,S2, ■■■, Sn}- 


N -1 

p(n,m-i,...,no) = En (si + i\G ni \si)(n i+ i\U Si \ni)p(so,no) (4) 


O- i=0 


The evolution operators for gene and protein states are: 
G ni = e-Wn.At, and Ugi = e -l H au with being the 

rate matrix for gene states transitions at fixed protein 
state ni while L Si is the Fokker-Planck operator or pro¬ 


tein diffusion corresponding to a fixed gene state S; for 
times ti . Explicitly the short time protein state propa¬ 
gator can be approximated by an Onsager-Machlup like 
Lagrangian function: 


(n i+ i|e Ls i Ati \ni) = J\T(Ati)exp 


1 / m+i-rii 

2 D Si (ni)\ A ti 



dni 


(5) 


The term A/"(AU) takes care of normalization and we have 
omitted the Jacobian and gauge terms in the Lagrangian 
that come from the nonlinearity of variable transforma¬ 
tions and the multiplicative nature of noise. These cor¬ 
rections are vanishingly small in the small noise limit 


which we are going to take late r 42,54 . The terms A Si and 
D Si stand for drift and diffusion coefficients correspond¬ 
ing to gene state s, at time U. We may formally take the 
limit of AU —» 0 and write the transition probabilities 
between protein numbers in a path integral form: 


N-l 


p(n,m- I,...,n 0 ) =E II ^ 




.A ti 


0 —Ati£ Si (n 


1 ) - 


2=0 


P dtC s (n.n,i) 


( 6 ) 


We consider the full trajectory probabilities rather 
than two time probabilities because after averaging 
over the gene states in Eq (H]) the reduced probability 


p{n, ni-i ,..., no) no longer obeys Markovian rules when 
the dichotomous noise has a finite correlation time. The 
eigenvalues A m (nj) of the gene state change rate matrix 
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W ni are Ao = 0 which corresponds to a steady state and 
Ai = — (k on + k Q ff) , the next largest eigenvalue which 
sets the time-scale Td = Aj -1 for gene state equilibration. 
The limit of Td —> 0 corresponds to the adiabatic regime 
while Td —> oo corresponds to the non-adiabatic regime of 
gene expression. To investigate how the timescale Td of 



(A) (B) 


FIG. 1. (A) Pictorial view of how dichotomous noise gener¬ 

ates churning cycles in the regime with finite gene switching 
Td (B) Top figure shows the dynamics with full noise (FN) 
generated by a master equation which accounts for diffusive 
motion of protein population and discrete jump transitions of 
the gene. In the bottom figure (DMN) the diffusive dynamics 
is replaced by deterministic equations of motion, while gene 
transitions are generated from the master equation while cou¬ 
pled to its deterministic counterpart. 


the dichotomous noise affects the trajectory probabilities 
we can perform a cumulant expansion of the exponential 
average (J6]). 


P dtC s (n,n.t) 

(3 J n 0 v ’ ' ' 


= exp 



dtC s (h, n, t) 


a 


k—oo 


+ E 



(7) 

We can see that in the adiabatic limit of fast gene switch¬ 
ing (Td —> 0) the higher order cumulants will vanish be¬ 
cause the state of the gene is instantaneously equilibrated 
at each time. In this regime the trajectory probability 
p(n, rii-i,no) can be approximated by the first cu¬ 
mulant term which corresponds to gene state averaged 
stochastic dynamics of protein numbers. Put another 
way the gene states cease to be correlated at subsequent 
time intervals At hence the probability p(n, rii-i, ...,no) 
acquires a markovian property and is therefore reduced 
to an effective Fokker-Planck description (Fig. |T] )■ 
The correlation function of the dichotomous gene vari¬ 
able — is given by Cx)Mjv(At) = ^exp(-^j), where 
Td = i— Xkjj the gene switching timescale and D = 

konkoff T d(9on — 9off) 2 is the strength of the dichoto¬ 
mous noise. The explicit expression for D in the context 
of stochastic gene expression was earlier obtained by^ 
where the term ’’churning noise” was coined, emphasiz¬ 
ing that cycles in the space of gene state and protein 
number variables generated by alternating seemingly fu¬ 
tile processes of binding/unbinding and protein produc¬ 
tion/degradation but that nevertheless these cycles leads 
to dispersion in the probability space. In the case of 
very slow gene switching (r^ —» oo) the memory of a 
gene state is preserved for a long time thus the trajec¬ 
tory probability can be expressed by a sum of two terms 


/ Jlwt , dtC s (h,n,t)\ \ dtCo . P’T ,d.tC\ 

(e "o (1 o) ) ps poe o( o) +pie Jn o (t t> ) cor¬ 

responding to diffusive motion of protein numbers with 
the initial state of the gene being ieither OFF or ON 
with probabilities po and pi respectively. In this extreme 
non-adiabatic regime one can simply model the dynamics 
of the gene circuit via Fokker Planck equations for each 
gene state independently. 

In the intermediate regimes (with finite Td), however, 
a significant component of a gene circuit’s stochasticity 
comes from cyclic motion in the protein number space 
and one needs to consider contributions from trajectories 
that perform such driven cyclic drift (Fig. |T| ). The 
contribution of cyclic motion is contained in the higher 
order cumulants. We now outline a particularly simple 
recipe for numerically sampling stochastic trajectories in 
the intermediate regime of gene switching which accounts 
for the cyclic motion in the protein number space. 

So far, other than the diffusive assumption for species 
represented by continuous population, no major approx¬ 
imations have been invoked in deriving the expression 
for the short time propagator for the coupled white- 
dichotomous noise dynamics of genetic circuit. We can 
further approximate of gene path sum by taking the limit 
of small birth/death noise D(n) —>• 0, which removes the 
protein fluctuations while preserving only the cyclic com¬ 
ponent of motion (Fig 1). In this limit the path inte¬ 
grals for each realization of gene state trajectory a in 
Eq[G]are replaced with deterministic trajectories and the 
whole path sum is dominated by the stochastic trajecto¬ 
ries which follow piecewise jump dynamics in the space 
of protein numbers: 


dn 

dt 

dp n (s,t) 


dt 


A(s, n) 


W nPn(s,t) 


(8a) 

(8b) 


Eq (SJr which governs protein dynamics completely ig¬ 
nores fluctuations in protein number that arise from 
birth/death processes. One may account for birth/death 
noise by adding a fluctuating linear term correction as 
noise B s i](t) that stems from the size expansion approx¬ 
imation^ to each deterministic mode A s . 


dn ., . „ . . 

— = A(s,n) +B s r](t ) 

(9a) 

dp it t} -«»?»<»•<> 

(9b) 

(r ] (t)r,(t')) = S(t-t') 

(9c) 


The system of equations [9] is one of the key results of 
this work focusing on gene networks. We note that piece- 
wise deterministic equations and other hybrid equations 
similar to Eq. [5] and Eq. [9] have already been used in 
the study of stochastic ion channels, neuronal dynamics, 
chemical reactions and finance 56,5 '. An interesting appli¬ 
cation of a hybrid scheme to model stochastic bursting 
effects in toggle switches can be found in the work by 
Bokes et al^S . Particularly attractive hybrid schemes for 
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studying general reaction systems can be found in the 
works of Salis-Kaznessis^ and Haseltine-Rawlings^S. In 
those schemes the master equation of the entire system is 
partitioned into fast and slow reaction blocks which are 
respectively modeled with a jump and continuous state 
Markov processes. Moreover, the partitioning of the re¬ 
actions is dynamically adjusted over the course of sim¬ 
ulations in order to achieve high quantitative accuracy. 
These and many other hybrid schemes while differing in 
the algorithmic details have the same underlying theme 
of exploiting time scale separation present in reactions 
for constructing physically and numerically sounds ap¬ 
proximations. The distinguishing features of the dichoto¬ 
mous noise based models proposed in this work are the 
emphasis on conceptual simplicity and the gene centric 
treatment of reactions where dichotomous noise is used 
exclusively for reactions involving gene states regardless 
of the rates of other reactions and white noise is used 
to model fluctuations in the molecular populations of all 
other species. Thus the proposed scheme is specifically 
aimed at genetic networks as a relatively coarse-grained 
(compared to kinetic Monte Carlo and other stochastic 
schemes) yet fairly simple explorative tool for studying 
steady state attractors and the associated intrinsic noise 
caused by the discrete dynamics of the promoter architec¬ 
ture (see the section below for the numerical algorithm). 
In the context of genetic networks we would like to single 
out the related works by Ge et alM who recently derived 
a fluctuation model for the self regulating gene model to 
study the kinetics of barrier crossing in the intermedi¬ 
ate regime of gene switching. Their fluctuation model is 
a particular application of Eq [8] It is also worth men¬ 
tioning the works by Karmakar-Bose^i and Zeiser et al^ 
who derived closed analytic expressions for moments and 
steady state distributions for hybrid models of self regu¬ 
lating gene. In the subsequent sections we will be refer¬ 
ring to the sets of equations [5] and [9] as the dichotomous 
Markov noise (DMN) and the dichotomous + linear noise 
(DMN+LNA) approximations respectively. Simulating 
DMN or DMN+LNA schemes can offer a computation¬ 
ally efficient alternative to direct kinetic MonteCarlo sim¬ 
ulations of the full master equation. Lastly the noise ob¬ 
tained via either the DMN or the DMN+LNA scheme can 
be used in order to decompose the total noise into genetic 
and non-genetic contributions via a 2 ot = <J 2 DMN + 
where cr 2 DMN is the variance that can be generated by the 
DMN scheme, a 2 ot is the variance computed at the level 
of master equation and & 2 BD is the non-genetic noise com¬ 
ing from birth-death type of reactive events. One may 
also find it useful to further decompose the birth death 
noise &bd =ij dmn+lna + a corr into a gaussian uncorre¬ 
lated part & BMN + LN a an d higher orders of noise <J 2 orr . 


III. HYBRID MULTI-SCALE STOCHASTIC 
SIMULATION OF GENETIC CIRCUITS 

The set of equations [8] and [9] offer a multi-scale scheme 
for simulating genetic networks, where one accounts for 
the dichotomous gene noise at the master equation level, 
while generating the birth death noise either via de¬ 
terministic dynamics in between stochastic gene state 
changes (DMN) or via a linear noise level of approx¬ 
imation (DMN+LNA). There are many ways of simu¬ 
lating a set of piece-wise deterministic equations and 
here we use a method similar to the one proposed in 
the work of Alfonsi et al^. The state of a gene cir¬ 
cuit (n, s) is represented by vectors of continuous (n) 
and discrete species (s). Since the deterministic protein 
evolution is coupled to stochastic gene switching, the rate 
coefficients (. k?j(n(t ))) of stochastic switching (s* —> Sj ) 
become time dependent and one can not use a conven¬ 
tional Gillespie type Monte Carlo algorithm. The numer¬ 
ical procedure for simulating hybrid stochastic schemes 
(DMN/DMN+LNA) outlined in the section II consists of 
the following steps. First, one specifies the initial condi¬ 
tions for protein and gene states (n 0 , s 0 ) for t = 0. Then 
the dynamics is evolved by solving either an ODE (Eq[51i 
) or SDE (Eq [SJi) sets with a reasonably small time step 
At which needs to be comparable or smaller than the typ¬ 
ical time for gene switching = (k on + fc 0 //) _1 - All the 
simulations in the paper are carried out with At = OTr^. 
Given the initial conditions and the specified time step 
the trajectories are generated by numerically integrat¬ 
ing the ODEs or SDEs for each time ti = i * At, with 
i = [0,1,2...]. The numerical integration continues until 
the condition for the stochastic gene switching (s* —> Sj) 
is reached. The switching condition is obtained by noting 
that the waiting time probability distribution for stochas¬ 
tic switching is given via P(r) = ea :p[ — J* +T a(n, t)dt) 
where a = Kj is the escape rate from the state 

(npSj) and r is the elapsed time after which stochastic 
event occurs. One can sample the stochastic transition 
events by generating a uniform random number r\ on the 
interval [0,1] and setting the condition for the stochastic 
switching as rq = exp[—f t T a(n, t)dt). If one has multi¬ 
ple genes in the system the criteria for choosing which one 
happens at the switching time is determined by drawing 
a second random uniform number r 2 and testing for the 
condition &tj < a-rq < Xqii Mj which determines 

the Si -A Sj gene switching event. 


A. Self regulating gene 

In this section we study the simplest model of the self 
regulating gene in order to explore and understand the 
ways white and dichotomous noise sources contribute to 
the overall noise level and also to test the quality of DMN 
and DMN+LNA approximations. The self-regulating 
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gene circuit involves the following reactions: 


ON OFF + P 

(10a) 

OFF + P ON 

(10b) 

P A- 0 

(10c) 

0 P 

(10d) 

0 END, p 

(lOe) 


The first two reactions represent the gene switching pro¬ 
cesses which are responsible for the dichotomous noise, 
while the other reactions correspond to birth and death 
events. The deterministic model for birth/death dy¬ 
namics is fi = g(s) — kn 7 where s stands for the gene 
state. The deterministic steady state attractors lead to 
differing protein numbers n = for each gene state 
s. The linear noise approximation around each attrac¬ 
tor is dn = f(s,n)dt + ^-dW where dW stands for 
the Wiener process. Varying the off rate k a ff while 
maintaining the relationship k on = 0.1 k a ff changes the 
adiabaticity with of the circuit smoothly going from 
the adiabatic to the non-adiabatic regime while not af¬ 
fecting the ’’equilibrium” properties of the binding pro¬ 
cess. To quantify the protein noise level we compute 
the Fano factor a 2 / g. The Fano factor quantifies the 
deviation from Poissonian stochastic kinetics which cor¬ 
responds to having a steady state distribution of protein 
numbers with a 2 /g = 1. The computed Fano factors 



FIG. 2. Dependence of Fano factor on the unbinding 
rate of promoter (fc 0 //) computed with dichotomous noise 
(DMN), DMN with linear noise (DMN+LNA) simulations 
and full noise simulations (FN) obtained by solving master 
equation. Inset shows the schematic topology of the model 
of self-regulated gene. Sample trajectories are shown on the 
left side. The values of the remaining rate coefficients are 
g 0 n = 25, g of f = 60 and km 1. 

shown in Fig. [2] indicate that both the DMN and the 
DMN+LNA schemes do quite well in capturing the order 
of magnitude of the noise in both the adiabatic and non- 
adiabatic limits. This is quite remarkable, since the noise 
in DMN comes entirely from the local steady nonfluctuat¬ 
ing states. This level of approximation does not account 
for the smallness of protein populations. Nevertheless 


the approximation provides a good estimation of the ob¬ 
servable Fano factors. The noise decomposition scheme 
outlined at the end of section II shows how noise from 
protein fluctuations caused by pure birth/death events 
compares to the total noise produced by the gene cir¬ 
cuit as a whole: cr% D /cr 2 ot ~ 16% in the non-adiabatic 
regime (fc Q // = 0.1s -1 ) and only cr 2 BD /a 2 ot ~ 70% even 
in the adiabatic regime (k a ff = 10s -1 ) for the present 
model circuit. Birth/death fluctuations become more im¬ 
portant as the circuit approaches the adiabatic regime. 
Nevertheless we see that the contribution of dichotomous 
noise remains significant throughout this wide range of 
kinetic regimes. 



FIG. 3. Comparison of probability distributions obtained 
with full master equation (red lines) with those of (A) DMN 
and (B) DMN+LNA schemes for different values of adiabatic¬ 
ity parameter k 0 ff. 

Inspection of the detailed features of steady state prob¬ 
ability distributions on the other hand reveals that pure 
dichotomous noise captures the skeleton of the distri¬ 
bution (Fig 3A), properly predicting whether the sys¬ 
tem’s distribution is bimodal or unimodal along with 
even the relative weight of attractor states. Neverthe¬ 
less birth/death induced fluctuations in each steady state 
attractor are needed to give the distribution more rea¬ 
sonable widths. The peak heights are overestimated in 
the DMN approximation since the system does not fluc¬ 
tuate much in protein number around the steady state 
attractors (go/k and gi/k). The DMN+LNA scheme on 
the other hand does a much better job at reproducing 
the distribution with higher accuracy. In the adiabatic 
limit both the DMN and the DMN+LNA schemes do 
quite well, although pure dichotomous noise alone nat¬ 
urally leads to more narrowly peaked distributions. A 
more stringent comparison between probability distribu¬ 
tions is provided by the Kullback-Leibler (KL) divergence 
Dkl(p\q ) = ^2 x p(x)log ^^ji also known as the relative 

entropy in information theory^. Dkl is a measure of 
statistical distinguishability between two probability dis¬ 
tributions p{x) and q(x). It is a non-symmetric and pos¬ 
itive function which vanishes only when the distributions 
are identical p{x) = q{ x). When p{x) is the exact distri¬ 
bution directly obtained from the master equation and 
q(x) is the approximate estimation via DMN + LN A or 
DMN schemes DKL{p\q) quantifies the amount of infor¬ 
mation in bits (see Eq. 11 which is irrevocably lost due 
to approximations. From the comparison (Fig. |3j) we see 
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FIG. 4. The Kullback-Leibler divergence Dkl(p\c[) be¬ 
tween exact distribution p and approximate estimations q 
for steady state probabilities obtained via (A) DMN and (B) 
DMN+LNA simulations plotted against the adiabaticity pa¬ 
rameter k 0 ff. 

that DMN+LNA is a better approximation in the inter¬ 
mediate range of gene switching and becomes worse in 
the very extreme limits. Thus the intermediate regime 
of gene switching is where such schemes are expected to 
shine. Many biochemical genetic circuits are expected to 
operate in this regime. 

Overall these results are encouraging in that they show 
that one can obtain the steady state landscapes of large 
scale gene networks by only including dichotomous gene 
noise while modeling birth death events with the simple 
additive white noise. Simulating every reaction at a sin¬ 
gle molecule level is not only computationally prohibitive 
but probably also unnecessary for most cases. Purely de¬ 
terministic models by themselves however are too crude 
for modeling gene networks. For instance the macro¬ 
scopic model of a self regulating gene with monomeric 
binding does not predict the possibility of bistable be¬ 
havior in the non-adiabatic regim o 45 ! 65 , hence including 
gene states explicitly is crucial for exploring steady states 
of genetic circuits. 

B. Extended toggle switch model and the importance of 
extinction states. 

In this section we consider a more complex and realistic 
circuit, the toggle switch. Toggle switches are ubiquitous 
bistable control modules. They have been implicated 


in regulating major cell fate decisions such as differen¬ 
tiation, cell death and development. Some well known 
examples of genetic toggle switches include A phage’s 
lysis/lysogeny switch, ccll-cycle control circuits, signal 
transduction and stem cell differentiation pathways, etc. 
Here we will look closely at an extended toggle switch 
circuit model proposed by Strasse r 66 : 67 , which displays 
interesting multi-attractor dynamics. The circuit of the 
extended toggle switch consists of two genes A and B. It 
is based on a two state gene expression model with ex¬ 
plicit transcription (producing mRNAs Ma and Mb) and 
translation (producing Pa and Pb) steps. The proteins 
Pa and Pb are mutually inhibiting which gives rise to a 
bistable behavior. The reaction scheme of the extended 
toggle switch circuit contains the following elementary 
steps: 


ON a 

OFF a + P B 

(11a) 

OFF a + Pb ON a 

(lib) 

ON b 

^4 OFFb + Pa 

(11c) 

OFF b + Pa ON b 

(lid) 

OFF a -^4 

■ OFFa + mRN A a 

(lie) 

OFF b ^ 

OFFb + mRN Ab 

(Ilf) 

mRN Aa 

^ mRN A A + P A 

(Hg) 

mRNAs 

mRN A b + P B 

(Ht) 


mRNA A ^ 0 

(in) 


mRNAs 0 

(nj) 


p ^4 0 

(ilk) 


P ^ 0 

(in) 

(llm) 


Reactions (a-d) are gene switching events, (e,f) corre¬ 
spond to transcription process, (g,h) to translation pro¬ 
cess and the remaining reactions describe the degrada¬ 
tion of mRNA and protein molecules. The deterministic 
equations describing the birth-death processes are: 


M a = fi(s A ,M A ) = s a /3a - 1 aM a (12a) 

Mb = f-2{sB,M B ) = sbPb - 'JbMb (12b) 

Pa = h(sB,MA, Pa) = QaMa - SaPa + k 0 ff B sb — k on , B (1 — sb)Pa ( 12 c ) 

Pb = U(sa,M b ,Pb ) = 9bM b - SbPb + k 0 ff A SA — k onA { 1 — sa)Pb (12d) 


We study the case of equal rate coefficients for genes A The variable s stands for the gene state of the first gene 
and B, making the toggle switch symmetric. (0 for OFF and 1 for ON ) and s' denotes the state of 

M S3 (s) = — (13a) 

1 

gM ss { S ') + k off (l- S ) 

+ k on • s 


P S8 (s) 


(13b) 














the second gene. The DMN+LNA equations are: 


J 


dM A (t) = fi(s A ,M A )dt + SAyffijidWi - \J-j A • M ss (s A )dW 2 (14a) 

dM B (t) = f 2 (s B ,M B )dt + ss^fpsdWz - \J^ B ■ M as (s B )dWi (14b) 

dP A (t ) = f 3 (s B ,M A , P A )dt + \Zg A M ss (s A )dW 5 - \JS A ■ P ss (s B )dW & (14c) 

dP B (t) = U(s A , M B ,P B )dt + sjg B M ss (s B )dW 7 - ^/S B ■ P ss (s A )dW 8 (14d) 


The stochastic transitions of the gene state s are again 
governed by a master equation. In the present model 



FIG. 5. Dependence of Fano factor on the mRNA 

degradation rate 7 computed with full noise (FN) dichoto¬ 
mous noise(DMN) and DMN with linear noise correction 
(DMN+LNA). Inset shows the schematic topology of the 
model of extended toggle switch. Sample trajectories are 
shown on the left side. The values of the remaining rate 
coefficients are g on = 0.05, <? 0 // = 0, S = 0.005, k on = 5, 
k 0 ff = 0.1 and /3 = 0.05 for both A and B genes. 

the basal rate of protein production in the repressed 
bound state has been set to zero ( g 0 ff = 0 ), which cre¬ 
ates strongly bistable behavior between extinction and 
steady production. The extinction attractors correspond 
to small or virtually non-existent mRNA and protein 
populations. Such extinction states are interesting to 
study with our methods both for understanding the how 
different types of noise sources impact the stability of 
bistable switches and also for finding the true limita¬ 
tions of our approximation schemes. Both the DMN and 
DMN+LNA equations are based on deterministic evolu¬ 
tion of molecular populations and they either completely 
neglect the birth/deatlr fluctuations or treat them ap¬ 
proximately. It is interesting to see how such schemes 
would behave in the regimes where the birth/death fluc¬ 
tuations are significant. By varying the mRNA degra¬ 
dation rate 7 we vary the average lifetime of mRNAs 
which controls the burst size of proteins. Thus we go 
from the regime of high birth/death noise (small 7 and 
~ 10 molecules) to low birth/death noise (large 7 and 
~ 10 2 molecules). The dependence of the Fano factor 
on 7 shows (Fig. [5]) that indeed both the DMN and 
DMN+LNA schemes capture the total noise levels of the 




FIG. 6 . Steady state probability distribution projected in the 
space of q A = N(mRNA a ) — N(mRNA b) and q B = N(P A )— 
N(Pb) order parameters reveals the multi-attractor nature of 
extended toggle switch. Comparing (A) Full Noisewith (B) 
DMN+LNA approximations for the different values of mRNA 
lifetime 7 . 


circuit quite well with the agreement improving steadily 
with the decreasing contribution of the birth/death noise. 
Detailed inspection of the probability distributions how¬ 
ever reveals that the DMN+LNA scheme fails to cap¬ 
ture metastable near-extinction states. To better visual¬ 
ize all of the steady state attractors we have projected 
the 2D probability distribution on the the order param¬ 
eters q A = M a — M b and q B = P A — P B which em¬ 
phasize the importance of states with M A « M B ~ 10 
and P A « P B ~ 10, where both proteins and mRNAs 
are present in equally low quantities. From the probabil¬ 
ity distribution at high birth/deatlr noise limit obtained 
from the simulation with full noise(see the top figure in 
Fig. [ 6 j\) we see that the toggle switch displays multi¬ 
attractor behavior with two main attractors correspond¬ 
ing to two configurations of the toggle switch with gene 
states s A = ON (OFF) and s B = OFF (ON) and the 
metastable attractor corresponding to both genes being 
simultaneously repressed (s^ = s B = ON) which gives 
rise to a narrow peak in the center of 2D probability dis¬ 
tribution. 

The steady state probability distributions obtained 
with the DMN+LNA scheme capture the main attrac- 
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tors in the system but do not reveal the presence of 
metastable states. This is because in the linear nose 
approximation of the birth death noise the fluctuations 
of discrete species are treated as uncorrelated brownian 
noise terms. Thus the correlated fluctuations that are 
obtained by treating the noise at the full master equa¬ 
tion level are essential to the existence of the extinction 
state attractors in this system. Such metastable attrac¬ 
tors while conceptually interesting have a negligible con¬ 
tribution to the total noise level. They are only observed 
when the proteins and mRNAs are present in extremely 
low quantities on the order of ~ 10. Once the numbers of 
molecules reach levels of ~ 100 the DMN+LNA scheme 
provides an excellent approximation to the full master 
equation as can be seen from the comparison of distribu¬ 
tions on the lower panel of Fig. [6] If one is interested in 
studying such metastable states with hybrid stochastic 
techniques, one possibility would be to treat the reac¬ 
tions involving extremely small numbers molecules into 
the space of discrete variables and treating them at the 
same master equation level as the gene states. 



50 300 550 800 1050 1300 1550 


N P 

FIG. 7. Residence times in the ON state of a gene vs the 
protein population size of an extended toggle switch model. 

Finally we study how the molecular population size af¬ 
fects the residences times of a gene being in the ON state 
which provides a measure of stability of an attractor. We 
vary the population size by changing the translation rate 
( g = gA = g B ) (Fig. [7]) for both genes. In the case of 
the full noise the residence times grow linearly with the 
population size (number of Proteins and mRNAs) which 
is in agreement with the previous studies on the extended 
toggle switch circui t 66 ! 67 . Overall due to the white noise 
contribution the noise at the DMN+LNA level allows 
for more fluctuations in the gene state which in turn re¬ 
sults in more transition events. This aspect allows faster 
exploration of gene circuits distribution, a feature that 
may be useful for quick scanning for the steady state at¬ 
tractors in large genetic networks. The residence times 
computed with the DMN scheme are obviously insensi¬ 
tive to the changes of protein numbers. This insensitivity 
is due to the deterministic treatment of these reactions 
involving species other than genes which do not affect the 


stochastic gene switching events. 

C. Conclusion 

In this work we have explored hybrid stochastic 
schemes based on treating the genetic noise of DNA bind¬ 
ing explicitly while making linear noise approximations 
for the birth and death reactions of other species. The 
theoretical motivation for such schemes along with the 
key assumptions and connections of hybrid schemes with 
the underlying master equation are elucidated in a path 
sum based analysis of trajectory statistics. We have ex¬ 
plored the steady state attractors of two classic toy mod¬ 
els of gene circuits, the self-regulating gene and an ex¬ 
tended version of a genetic toggle switch. Our findings 
indicate that gene state change noise accounts for a large 
fraction of the total noise for a wide range of time-scales. 
The noise decomposition obtained by comparing the total 
noise estimated by the different schemes for the extended 
model of toggle switch circuit revealed that the noise in¬ 
duced metastable attractor originates from correlations 
of discrete molecular populations. Our study suggests 
that stochastic hybrid schemes will be useful tools for ex¬ 
ploring large scale genetic networks, which display com¬ 
plex multi-attractor dynamics caused by discrete fluc¬ 
tuations in the promoter architecture and should find 
their place alongside more traditional full scale kinetic 
Monte Carlo simulations. For future work it would be 
interesting to use such hybrid schemes to model the ef¬ 
fects of gene co-localization, heterogeneous distribution 
of transcription factors in nuclear/cytoplasmic milieu and 
1D/3D diffusional motions by adding additional deter¬ 
ministic steps for each molecular component and enlarg¬ 
ing the number of internal discrete promoter states. The 
latest studies have shown that in eukaryotic cells gene 
expression is an intricate function of chromatin config¬ 
uration as well as spatial distribution of decoy binding 
sties on the genome. Thus it would be exciting to use 
hybrid schemes accounting for internal gene states with 
a high level of detail via a discrete web of Markov states 
to get a good understanding of some of these complex 
facets in eukaryotic gene regulation. At last a possible 
extension to hybrid schemes with coupled diffusion intro¬ 
duced for selected species could also be an attractive way 
of treating in a coarse grained manner the cases where 
the well stirred assumption does not hold. 
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